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Improving convergence in smoothed particle hydrodynamics 
simulations without pairing instability 
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ABSTRACT 

The numerical convergence of smoothed particle hydrodynamics (SPH) can be severely re- 
stricted by random force errors induced by particle disorder, especially in shear flows, which 
are ubiquitous in astrophysics. The increase in the number Nh of neighbours when switching 
to more extended smoothing kernels at fixed resolution (using an appropriate definition for 
the SPH resolution scale) is insufficient to combat these errors. Consequently, trading reso- 
lution for better convergence is necessary, but for traditional smoothing kernels this option is 
limited by the pairing (or clumping) instability. Therefore, we investigate the suitability of the 
Wendland functions as smoothing kernels and compare them with the traditional B-splines. 
Linear stability analysis in three dimensions and test simulations demonstrate that the Wend- 
land kernels avoid the pairing instability for all Nh, despite having vanishing derivative at the 
origin (disproving traditional ideas about the origin of this instability; instead, we uncover a 
relation with the kernel Fourier transform and give an explanation in terms of the SPH density 
estimator). The Wendland kernels are computationally more convenient than the higher-order 
B-splines, allowing large Nh and hence better numerical convergence (note that computational 
costs rise sub-linear with Nh)- Our analysis also shows that at low Nh the quartic spline kernel 
with Nh ~ 60 obtains much better convergence then the standard cubic spline. 

Key words: hydrodynamics — methods: numerical — methods: A^-body simulations 



1 INTRODUCTION 

Smoothed particle hydrodyna mics (SPH) is a particl e -based nu- 
merical met hod, pioneered by iGingold & MonaghanI d 19771) and 
lLucvl ( ll977h . fo r solving the equations of hydrodynami c s (recent 
revie^ys include iMonaghai] |2005| , 1201 2l : iRosswogf [20091 : ISpringell 
l201(]| : |Pricel2012h . In SPH, the particles trace the flow and serve as 
interpolation points for their neighbours. This Lagrangian nature of 
SPH makes the method particularly useful for astrophysics, where 
typically open boundaries appl y, though it beco mes increasingly 
popular also in engineering (e.g. lMonaghanll20 1 2h . 

The core of SPH is the density estimator: the fluid density is 
estimated from the masses m, and positions x, of the particles via 
(the symbol ' denotes an SPH estimate) 



tion fo r the particles are derived, following iNelson & Papaloizoul 
( Il994l) . via a variational principle from the discretised Lagrangian 

£ = i:,m,[^x':-u(f,„sd] (2) 

jMonaghan & Pried lioOlh . Here, u(p,s) is the internal energy as 
function of density and entropy s (and possibly other gas proper- 
ties), the precise functional form of which depends on the assumed 
equation of state. The Euler-Lagrange equations then yield 



nij dxi 



-^V.W{x,j,h,) + -^V,W(x,j,hj) 



where Xjj = x, - Xj and P, = du/dpi, while the factors 



p(Xi) ~Pi = W(Xi-Xj,hj), 



(1) 



1 dihip,) 



1 



(3) 



(4) 



where W{x,h) is the smoothing kernel and hj the smoothing scale, 
which is adapted for each particle such that /ijp, =constant (with 
V the number of spatial dimensions). Similar estimates for the 
value of any field can be obtained, enabling discretisation of the 
fluid equations. Instead, in conservative SPH, the equations of mo- 
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vKpj din hj 

dSpringel & Hemauisdl2002l:lMonaghanll2002h arise from the adap- 
tion of hi ( Nelson & PapaloizoiJ) such"that /jj'p, =constant. 

Equation ([3j is a discretisation of px = -VP, and, because 
of its derivation from a variational principle, conserves mass, lin- 
ear and angular momentum, energy, entropy, and (approximately) 
circularity. However, its derivation from the Lagrangian is only 
valid if all fluid variables are smoothly variable. To ensure this, 
in particular for velocity and entropy, artificial dissipation terms 
have to be added to x, and ;(,. Recent progress in restricting 
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such dissipation to regions of compressive flow dCuIlen & DehnenI 
l2010l : lRead & HavfieldlioT^) have greatly improved the ability to 
model contact discontinuities and their instabilities as well as near- 
inviscid flows. 

SPH is not a Monte-Carlo method, since the particles are not 
randomly distributed, but typically follow a semi-regular glass-like 
distribution. Therefore, the density (and pressure) error is much 
smaller than the > 15% expected from Poisson noise for ~ 40 neigh- 
bours and SPH obtains OQr) convergence. However, some level 
of particle disorder cannot be prevented, in particular in shear- 
ing flows (as in turbulence), where the particles are constantly re- 
arranged (even in the absence of any forces), but also after a shock, 
where an initially isotropic particle distribution is squashed along 
one direction to become anisotropic. In such situations, the SPH 
force ( [3) in addition to the pressure gra dient contains a random 'Eq 
error (Read. Havfield & Agert3l2010lR and SPH converges more 
slowly than O(li^). Since shocks and shear flows are common in 
star- and galaxy-formation, the 'Eq errors' may easily dominate the 
overall performance of astrophysical simulations. 

One c an dodge the 'Eq error' by using other discretisations of 
px - -VP ( lMorrisiri996l : lAbelll201ll) . However, such approaches 
unavoidably abandon momentum conserv ation and hen ce fail in 
practice, in particular, for strong shocks (iMorrisI [19961) . Further- 
more, with such modifications SPH no longer maintains particle or- 
der, which it otherwise automatically achie ves. Thus, t he 'Eq error' 
is SPH's attempt to resunect particle order (Price I2OI2I) and prevent 
shot noise from afi'ecting the density and pressure estimates. 

Another possibility to reduce the 'Eo eiTor' is to subtract an av- 
erage pressure from each particle's P, in equation ([Sj- Efi'ectively, 
this amounts to adding a negative pressure term, which can cause 
the tensile instability (see 93.1.2b . Moreover, this trick is only use- 
ful in situations with little pressure var iations, perhaps i n simula- 
tions of near-incompressible flows (e.g. lMonaghanll20lTI) . 

The only remaining option for reducing the 'Eo error' appears 
an increase of the number of particles contributing to the den- 
sity and force estimates (contrary to naive expectation, the compu- 
tational costs grow sub-linear with N^). The traditional way to try 
to do this is by switching to a smoother and more extende d ker- 
nel, e nabling larger Nh at the same smoothing scale h (e.g. IPricd 
|20l3)- 

However, the degree to which this approach can reduce the 
'Eo errors' is limited and often insufliicient, even with an infinitely 
extended kernel, such as the Gaussian. Therefore, one must also 
consider 'stretching' the smoothing kernel by increasing h. This 
inevitably reduces the resolution, but that is still much better than 
obtaining erroneous results. Of course, the best balance between 
reducing the 'Eo error' and resolution should be guided by results 
for relevant test problems and by convergence studies. 

Unfortunately, at large N/j the standard SPH smoothing 
kernels become unstable to the pairing (or clumping) instabil- 
ity (a cousin of the tensile instability), when particles form 
close pairs reducing the efi^ective ne ighbour number. The pair- 
ing instability (first mentioned by ISchilBler & Schmit3 Il98lh 
has traditionally been attributed to the diminution of the re- 
pulsive force between close neighbours approach i ng each other 
I Schiifiler & Schrnit?. 'Thomas & Couchman' '1992', 'Herand ll994 
ISwegle, Hi cks & Attawav 1995, Springel 2010, Price 2o{l)- Such 
a diminishing near-neighbour force occurs for all kernels with an 
inflection point, a necessary property of continuously difi'erentiable 



kernels. Kernels without that property ha ve been proposed and 
shown to be more stable (e.g. iRead et all) . However, we provide 
demonstrably stable kernels with inflection point, disproving these 
ideas. 

Instead, our linear stability analysis in Section |3] shows that 
non-negativity of the kernel Fourier transform is a necessary con- 
dition for stability against pairing. Based on this insight we propose 
in Section |2] kernel functions, which we demonstrate in Section |4] 
to be indeed stable against pairing for all neighbour numbers N/j, 
and which possess all other desirable properties. We also present 
some further test simulations in Section |4] before we discuss and 
summarise our findings in Sections |5]and|6] respectively. 



2 SMOOTHING MATTERS 
2.1 Smoothing scale 

SPH smoothing kernels are usually isotropic and can be written as 

W{x,h) = h-"w(\x\/h) (5) 

with a dimensionless function w(r), which specifies the functional 
form and satisfies the normalisation 1 = Jd''jt:vv(|jc|). The re-scaling 
h — > ah and ivir) — » a'wiar) with a>0 leaves the functional form 
of W{x) unchanged but alters the meaning of h. In order to avoid 
this ambiguity, a definition of the smoothing scale in terms of the 
kernel, i.e. via a functional h = h[W{x)], must be specified. 

In this study we use two scales, the smoothing scale h, defined 
below, and the kernel -support radius H, the largest \x\ for which 
W(x} > 0. For computational efficiency, smoothing kernels used in 
practice have compact support and hence finite H. For such kernels 



W{x,h) = H-"w{\x\/H), 



(6) 



where w(r) = for r > 1 and w(r) > for r< I. H is related to the 
average number Nh of neighbours within the smoothing sphere by 



(7) 



with Vy the volume of the unit sphere. H and Nh are useful quan- 
tities in terms of kernel computation and neighbour search, but not 
good measures for the smoothing scale h. Unfortunately, there is 
some confusion in the SPH literature between H and h, either be- 
ing denoted by 'Ii and refeiTed to as 'smoothing length'. Moreover, 
an appropriate definition of h in terms of the smoothing kernel is 
lacking. Possible definitions include the kernel standard deviation 

(r- = v-' j'&xx^W{x,h), (8) 

the radius of the inflection point (maximum of IVWI), or the ratio 
M7|VW| at the inflection point. For the Gaussian kernel 



W{x) = N{0,a-) = 



1 



exp 



X- 



(9) 



all these give the same result independent of dimensionality, but not 
for other kernels ('triangular' kernels have no inflection point). Be- 
cause the standard deviation ([§) is directly related to the numerical 
resolution of sound waves 03.1.3t . we set 



h = la. 



(10) 



' Strictly speaking, the 'Eq enor' term of iRead et al.l is only the dominant 
contribution to the force errors induced by particle discreteness. 



In practice (and in the remainder of our paper), the neighbour 
number is often used as a convenient parameter, even though it 
holds little meaning by itself. A more meaningful quantity in terms 
of resolution is the average number A';, of particles within distance 
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Wendland V = 2, 3 i/rj . i = ( 1 - r)l ( 1 + 4r) 

Wendland , v = 2, 3 (/'4.2 = ( 1 - '")+ ( 1 + 6r + ^ ) 

Wendland C^ v = 2, 3 1A5.3 = (1 - r)^ (1 + 8r + 25)- + 32r3 ) 

Table 1. Functional forms and various quantities for the B-splines (equation II 11 and Wendland functions (equation II 2t in v = 1-3 spatial 
dimensions. (■)+ = max{0,-|. C is the normalisation constant, tr the standard deviation (equation |8], and h = 2(T the smoothing scale. Note that 
the Wendland functions of given differentiabihty are identic al for v = 2 a nd v = 3 but dift'er from those for v = 1 . 1^2,1 (the Wendland function 
in ID) has already been used in the second SPH paper ever jLucvll977l) . but for 3D simulations, when it is not a Wendland function. 



7 21_ _5_ _1_ 

n Ik 72 15 

9 495 7 2 

78 B65 J_ J_ 

7i 64;r 70 24 



— 1.897367 1.936492 

— 2.171239 2.207940 

— 2.415230 2.449490 



h, given by Ni, = (h/Hy'Nu for kernels with compact support, or the 
ratio hip/m)''" between h and the average particle separation. 



2.2 Smoothing kernels 

After these definitio ns, let us list the desirable prop erties of the 
smoothing kernel (cf. lFulk & Ouinnlll996l : IPricel2012h . 

(i) equation l|T) obtains an accurate density estimate; 

(ii) W(x,h) is twice continuously differentiable; 

(iii) SPH is stable against pairing at the desired N^', 

(iv) W(x,h) and VW(x,h) are computationally inexpensive. 

Here, condition Q implies that W{x, h) — > S(x) as /i — > but also 
that W(x,h) > is monotonically declining with \x\; condition ^ 
guarantees smooth forces, but also implies Viy(O) = 0. 



2.2.1 B-splines 



The m ost used SPH kernel functions are the ISchoenberi 
1946h B-spline functions, generated as ID Fourier transform: 



( Monaghan & Lattanziolll985h 



w(r) = C6„(r), b„(r) = 



sin k/n 
kin 



cos kr dk. 



(11) 



with normalisation constant C. These kernels consist of 
piece-wise polynomials of degree n-\ (see Table [T} and are n — 2 
times continuously differentiable. Thus, the cubic spline (n = 4) is 
the first useful, but the quartic and quintic have also been used. For 
large n, the B-splines approach the Gaussian: b„ — > N(0,H-/3n) 
(this follows from footnote |2]and the central limit theorem). 

Following Monaghan & Lattanzio, h = 2H/n is conventionally 
used as smoothing scale for the B-splines independent of v. This is 
motivated by their original purpose to interpolate equidistant one- 
dimensional data with spacing h, but cannot be expressed via a 



^ By this definition they are the ;i-fold convolution (in one dimensi on) of 
bi{r) with i t self (m odulo a scahng), and hence are identical to the llrwini 
)l927hjHaii )l927h probability density for the sum r of n independent ran- 
dom variables, each unifonnly distributed between -l/;i and \/n. 



functional h = h[W(x)]. Moreover, the resulting ratios between h 
for the b„ do not match any of the definitions discussed abov^ 

Instead, we use the more appropriate h = 2cr also for the B- 
spline kernels, giving H x 1.826/; for the cubic spline in 3D, close 
to the conventional H = 2h (see Table[T](. 



2.2.2 'Triangular' kernels 

At low order n the B-splines are only stable against pairing for mod- 
est values of Nu (we will be more precise in Section |3j, while at 
higher n they are computationally increasingly complex. 

Therefore, alternative kernel functions which are stable for 
large Nh are desirable. As the pairing instability has traditionally 
been associated with the presence of an inflection point (minimum 
of w'), functions w(r) without inflection point have been proposed. 
These have a triangular shape at r ~ and necessarily violate point 
(ii) of our list, but avoid the pairing instabilit For comparison we 
consider one of them, the 'HOCT4' kernel of lkead et al] ( l2010l) . 



2.2.3 Wendland functions 

The linear stability analysis of the SPH algorithm, presented in the 
next Section, shows that a necessary condition for stability against 
pairing is the non-negativity of the multi-dimensional Fourier trans- 
form of the kernel. The Gaussian has non-negative Fourier trans- 
form for any dimensionality and hence would give an ideal kernel 
were it not for its infinite support and computational costs. 

Therefore, we look for kernel functions of compact support 
which have non-negative Fourier transform in v dimensions and are 



' Fig. 2 of l Pried i2012l) seems to suggest that with this scaling the B-splines 
approach the Gaussian with it = Ii/^. However, this is just a coincidence 
f or » = 6 (quintic sp line) since cr = ^[nj\2h for the B-splines in ID. 
* [Xh omas & Couc hman ( 1992) proposed such kernels only for the force 
equation ^3)> t"Jl to keep a smooth kemel for the density estimate. However, 
such an approach cannot be derived from a Lagr angian and h ence necessar- 
ily violates energy and/or entropy conservation ( |Pricel2012b . 
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Gaussian 
H0CT4 

cubic spline 
quartic spline 
quintic spline 
Wendland C= 
Wendland C 
Wendland C= 




Figure 1. Kernels of Table[T] the Gaussian, and the HOCT4 kernel of Read 
et al. (2010) scaled to a common resolution scale of h = 2<t for 3D {top: 
linear plot, aiTows indicating \x\ = H\ bottom: logarithmic plot). 



low-o rder polynomial jj in r. This is precisely the defining property 
of the lWendlandl ( Il995h functions, which are given by 



with (■)+ = max{0, ■) and the linear operator 
I[f](r) = fsf(s)ds. 



(12) 



(13) 



In V spatial dimensions, the functions i//(i,(\x\) with C = k+l + Lv/2J 
have positive Fourier transform and are 2k times continuously dif- 
ferentiable. In fact, they are the unique polynomials in x\ of mini- 
mal degree with these properties (IWendlandlll995Ll2005h . For large 
k, they approach the Gaussian, which is the only non-trivial eigen- 
function of the operator I. We list the first few Wendland functions 
for one, two, and three dimensions in Table [T] and plot them for 
v = 3 in Fig.[T] 



2.3 Kernel comparison 

Fig-dplots the kernel functions w(r) of Table [T] the Gaussian, and 
the HOCT4 kernel, all scaled to the same h = 2a for y = 3. Amongst 
the various scalings (ratios for hjH) discussed in 32. 11 above, this 
gives by far the best match between the kernels. The B-splines and 
Wendland functions approach the Gaussian with increasing order. 
The most obvious difference between them in this scaling is their 
central value. The B-splines, in particular of lower order, put less 
emphasis on small r than the Wendland functions or the Gaussian. 

Obviously, the HOCT4 kernel, which has no inflection point, 
differs significantly from all the others and puts even more empha- 
sis on the centre than the Gaussian (for this kernel a x 0.228343/f). 



Gaussian 
H0CT4 
cubic spline 
quartic spline 
quintic spline 
Wendland 
Wendland 
Wendland C 



<t 0.001 T 




0.0001 r 



Figure 2. Fourier transforms W(k) for the Gaussian, the HOCT4 and the 
kernels of Table[T]scaled to the same common scale h = 2tT. Negative values 
are plotted with broken curves. 



2.4 Kernel Fourier transforms 

For spherical kernels of the form (|6), their Fourier transform only 
depends on the product H\k\, i.e. W(k) = w{H\k\). In 3D (!Fk denotes 
the Fourier transform in v dimensions) 

w{k) = T3 [w{r)] (k) = 4nK'' J^°° sin(;^f-) w(r) rdr (14) 

which is an even function and (up to a normalisation constant) 
equals -K^'dfi[w]/dK. For the B-splines, which are defined via 
their ID Fourier transform in equation (II 1) . this gives immediately 

[b„(r)] (k) = 3 )"^' sin" |(l - ^ cot ^) (15) 

(which includes the normalisation constant), while for the 3D 
Wendland kernels 

(1 d \**' 

(we abstain from giving individual functional forms). 

All these are plotted in Fig.|2]after scaling them to a common 
h = 2a. Notably, all the B-spline kernels obtain W <Q and oscil- 
late about zero for large k (which can also be verified directly from 
equationfTSt, whereas the Wendland kernels have W(k) > at all k, 
as does the H0CT4 kernel. As non-negativity of the Fourier trans- 
form is necessary (but not sufficient) for stability against pairing at 
large Nh (see 33.1.2t , in 3D the B-splines (of any order) fall prey to 
this instability for sufficiently large N^, while, based solely on their 
Fourier transforms, the Wendland and H0CT4 kernels may well be 
stable for all neighbour numbers. 

At large |^| (small scales), the HOCT kernel has most power, 
caused by its central spike, while the other kernels have ever less 
small-scale power with increasing order, becoming ever smoother 
and approaching the Gaussian, which has least small-scale power. 

The scaling to a common h = 2a'm Fig.[2]has the effect that the 
W{k) all overlap at small wave numbers, since their Taylor series 



(16) 



W(k)-- 



l-la^k' 



+ 0(\kf). 



(17) 



' Polynomials in would avoid the computation of a square root. How- 
ever, it appears that such functions cannot possibly have non-negative 
Fourier transfomi (H. Wendland, private communication). 



2.5 Density estimation and correction 

The SPH force Q is inseparably related, owing to its derivation 
via a variational principle, to the derivative of the density estimate. 
Another important role of the SPH density estimator is to obtain 
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Figure 3. SPH density estimate (T) obtained from particles in three- 
dimensional densest-sphere packing {solid), glass (squares), or pairing 
(crosses), plotted against neighbour number Nh for the kernels of Table [T] 
and the HOCT4 kernel (colour coding as in Figs.[T]§J3). For the Wendland 
kernels, the corrected density estimate ( 1181 is shown with dashed curves 
and triangles for densest-sphere packing, and a glass, respectively. 



accurate values for P, in equation lO, and we will now assess the 
performance of the various kernels in this latter respect. 

In Fig. [3] we plot the estimated density (TJ vs. neighbour 
number Nh for the kernels of Table [T] and particles distributed in 
three-dimensional densest-sphere packing (solid curves) or a glass 
(squares). While the standard cubic spline kernel under-estimates 
the density (only values A'^ ;£ 55 are accessible for this kernel ow- 
ing to the pairing instability), the Wendland kernels (and Gaussian, 
not shown) tend to over-estimate it. 

It is worthwhile to ponder about the origin of this density 
over-estimation. If the particles were randomly rather than semi- 
regularly distribu ted, p obtained fo r an unoccupied position would 
be unbiased (e.g. ISilverman|[l986l) . while at a particle position the 
self contribution /)j,W(0,/7,) to p, results in an over-estimate. Of 
course, in SPH and in Fig.[3]particles are not randomly distributed, 
but at small the self-contribution still induces some bias, as ev- 
ident from the over-estimation for all kernels at very small N^- 

The H0CT4 kernel of Read et al. (2010, orange) with its cen- 
tral spike (cf. Fig.[TJ shows by far the worst performance. However, 
this is not a peculiarity of the HOCT4 kernel, but a generic property 
of all kernels without inflection point. 

These considerations suggest the corrected density estimate 



--pi-emiW(Q,hi), 



(18) 



which is simply the original estimate ([T) with a fraction e of the 
self-contribution subtracted. The equations of motion obtained by 
replacing p, in the Lagrangian (|2j with p,,corr are otherwise identical 
to equations l|3j and ^ (note that dih^p-^^^y/dlnhi = d{h]p^ld\nh,, 
since Kpi and /!)'p,_corr differ only by a constant), in particular the 
conservation properties are unaffected. From the data of Fig. [3] we 
find that good results are obtained by a simple power-law 



e = 6100 (A'/// 100)-" 



(19) 



with constants 6ioo and a depending on the kernel. We use (eioo, ce) 
= (0.0294,0.977), (0.01342,1.579), and (0.0116,2.236), respec- 
tively, for the Wendland C^, C"*, and kernels in v = 3 dimensions. 

The dashed curves and triangles in Fig.|3]demonstrate that this 
approach obtains accurate density and hence pressure estimates. 



3 LINEAR STABILITY AND SOUND WAVES 

The SPH linear stability analysis considers a plane-wave perturba- 
tion to an equilibrium configuration, i.e. the positions are perturbed 
according to 



Xi Xi + a exp {i[k-Xi - wf]) 



(20) 



with displacement amplitude a, wave vector k, and angular fre- 
quency (jj. Equating the forces generated by the perturbation to lin- 
ear order in a to the acceleration of the perturbation yields a disper- 
sion relation of the form 



aP(k) = clTu. 



(21) 



This is an eigenvalue problem for the matrix P with eigenvector a 
and eigenvalue oj^. The exact (non-SPH) dispersion relation (with 
c* = dPjdp, P = p'du/dp at constant entropy) 

c~ akk = oTa 



(22) 



has only one non-zero eigenvalue a) = c k with eigenvector a \\ k, 
con'esponding to longitudinal sound waves propagating at speed c. 

The actual matrix P in equation | |21| | depends on the details of 
the SPH algorithm. Fo r conservative SPH with equation of motion 
ll3t, lMonaghMil ( |2005i) gives it for P oc p^" in one spatial dimension. 
We derive it in appendix |A] for a general equation of state and any 
number v of spatial dimensions: 



2P 
P 



1 

2vp-Cl^ dlnli I 



(23) 



where x'^' is the outer product of a vector with itself, bars denote 
SPH estimates for the unperturbed equilibrium, t = pQu, and 



cos kXj)V^^^W{Xj,h), 



(24a) 



(24b) 



(24c) 



u(k) = —L- msmk-XiVW(Xi,h), 
p{n] ^ 

U(k)= mil- 

^_ 1 d\n(pCl) 
V d\nh 

Here and in the remainder of this section, curly brackets indicate 
terms not present in the cas e of a constan t /; = h, when our resul ts 
reduce to relations given bv lMon-isl(ll996l) and [Read et al.l(l2010h . 

Since P is real and symmetric, its eigenvalues are real and 
its eigenvectors mutually orthogonajf). The SPH dispersion relation 
Olb can deviate from the true relation ( I22t in mainly two ways. 
First, the longitudinal eigenvalue co^^ (with eigenvector ay || k) may 
deviate fro m c^k~ (wrong sound speed) or e ven be negative (pairing 
instability; lMorrislll996riMonaghanll2000l) . Second, the other two 
eigenvalues oj^, ^ may be significantly non-zero (transverse insta- 
bility for oj^i 2 < or transverse sound waves for oj^j ^ > 0). 

The matrix P in equation ( 123b is not accessible to simple in- 
terpretation. We will compute its eigenvalues for the various SPH 
kernels in § 33.21 -3 and Figs. l4l6( but first consider the limiting cases 
of the dispersion relation, allowing some analytic insight. 



3.1 Limiting cases 

There are three spatial scales: the wavelength A = 2n/\k\, the 
smoothing scale h, and the nearest neighbour distance £/„„. We will 



' If in equation (5) one omits the factors £1, but still adapts /j, to obtain 
Kpi =constant, as some practitioners do, then the resulting dispersion rela- 
tion has an asymmetric matrix P with potentially complex eigenvalues. 
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Figure 4. Stability of conservative SPH witli tlie Gaussian (truncated at 16it) and cubic spline kernel for densest-sphere packing and Pxp^^^: contours 
of (iij/c^t^ (top) and uj^^^lc^k- {bottom), the longitudinal and smallest transverse eigenvalues of P/c-k^ respectively, over wave number \k\ and smoothing 
scale h = 2cr (and Nh or Ni,) both scaled to the nearest-neighbour distance d„n. The left and right sub-panels are for t cc (1, 1, 1) (perpendicular to hexagonal 
planes), and fcoc(l, 1,0) (nearest-neighbour direction), respectively (other wave vectors give similar results). Red contours are for ai- <0 (implying the paiiing 
instability in the top panels) and are logaiithmically spaced by 0.25 dex. Blue contours are also logarithinically spaced between 10"* and 0.1, black contours 
are linearly spaced by 0.1, while good values for cj^ are cj^/c^k- = 0.95, 0.99, 0.995, 0.999, 1.001 1.005, 1.01, 1.05 (green), and 1 (cyan). The dashed line 
indicates a sound wave with wavelength A = 8// = 16o". For the Gaussian kernel |(i<^,| in the bottom panels is often smaller than our numeiical precision. 



separately consider the limit A»h of well resolved waves, the con- 
tinuum limit h » rf„„ of large neighbour numbers, and finally the 
combined limit A»h» d„„. 



3.1.1 Resolved waves 

If \k\h <g: 1, the argument of the trigonometric functions in equations 
(I24ti.b) is always small and we can Taylor expand thenj^. If we also 
assume a locally isotropic particle distribution, this gives to lowest 
order in \k\ (I is the unit matrix; see also 3A3t 



2v 



llU*'' + 



v+2 



with the eigenvalues 



^11 - g.2 I 5 ( 3v 



(25) 

(26a) 
(26b) 



The error of these relations is mostly dictated by the quality of the 
density estimate, either directly via p, c, and P, or indirectly via H. 
The density correction method of equation ( 118) can only help with 
the former, but not the latter. The difference between constant and 
adapted h is a factor 4/9 (for 3D) in favour of the latter. 



3.1.2 Continuum limit 

For large neighbour numbers N^, H'» d^„, ilScti — > I, p p and 
the sums in equations I l24t t.b) can be approximated by integr al^l 



u^-k W(k), and U ^ A:'^' 'W(k) 



(27) 



with W(k) the Fourier transform of W(x,h). Since W(k) = w(H\k\), 
we have dW/d\nh = fc-Vj,VK and thus from equation ( 123) 



^ IP I 
H/+— 1- 

pc^ \ 



W- 



I) 



(28) 



M'(O) = 1, but towards larger \k\ the Fourier transform decays, 
W{k) < 1, and in the limit h\k\ » 1 or /i « W ^ 0: short sound 
waves are not resolved. 

Negative eigenvalues of P in equation ( 128) . and hence linear 
instability, occur only if W itself or the expression within square 
brackets are negative. Since < 1, the latter can only happen 
if P < 0, which does usually not arise in fluid simulations (un- 
less, possibly, one subtracts an average pressure), but possibly in 
elasticity simulations of solids (Gray, Monaghan & Swift 20()J), 
when it causes the tensile instabi lity (an equivalent eff^ect is present 
in sm o othed-partic l e MH D, see lPhillips & Monaghanll 19851 ; IPricel 
l2012h . iMonaghaiJ ( |2000|) proposed an artificial repulsive short- 
range force, effectuating an additional pressure, to suppress the ten- 
sile instability. 

The pairing instability, on the other hand, is caused by W < 
for some H\k\ > kq. This instability can be avoided by choosing the 
neighbour number Nh small enough for the critical wave number kq 



^ In his analysis of ID SPH. lRasIcI )200C| ) also considers this simplification, 
but interprets it incoiTectly as the limit |fc|rf„n « 1 regardless of h. 



Assuming a uniform particle distribution. A better approximation, which 
does not require rfnn is to assume some radial distribution fimction 

g(r) (as in statistical mechanics of glasses) for the probability of any two 
particles having distance r. Such a treatment may well be useful in the con- 
text of SPH, but it is beyond the scope of our study. 
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Jfx(l 1,1) quartic spline k«(l,l,0) 



Jc«(l,l,l) Wendland Jc«(l,l,0) 




Figure 5. As Fig. |4] but for the quartic and quintic spline and the H0CT4 kernel (left) and the Wendland to kernels (right). The fine details of these 
contours are specific to the densest-sphere packing and will be different for more realistic glass-like particle distributions. 
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to remain unsampled, i.e. > K^ycp^st or H < d„„Ko/n (thiougti such 
small H is no longer consistent with the continuum limit). 

However, if the Fourier transform of the kernel is non-negative 
everywhere, the pairing instability cannot occur for large Nh- As 
pairing is typically a problem for large Nh, this suggests that ker- 
nels with W(k) > for every k are stable against pairing for all 
values of Nh, which is indeed supported by our results in 34.11 

3.1.3 Resolved waves in the continuum limit 

The combined limit of /I » /? » dnn is obtained by inserting the 
Taylor expansion ( I17t of W into equation ( I28t . giving 

d] = c^k^[\+crk^y-^[{2lv] + 1 - r] +0{h*\kt)) . (29) 

iMonaghanI j2005l) gave an equivalent relation for y = 1 when the 
expression in square brackets becomes 3 - y or 1 - y (for adapted or 
constant h, respectively), which, he argues, bracket all physically 
reasonable values. However, in 3D the value for adaptive SPH be- 
comes 5/3-7, i-S- vanishes for the most commonly used adiabatic 
index. 

In general, however, the relative error in the frequency is 
oc cp-k~ Bc (cr/A)-. This shows that h = 2cr is indeed directly pro- 
portional to the resolution scale, at least concerning sound waves. 



3.2 Linear stability of SPH kernels 

We have evaluated the eigenvalues ai^ and w^, j of the matrix P in 
equation i23i for all kernels of Table [T] as well as the H0CT4 and 
Gaussian kernels, for unperturbed positions from densest- sphere 
packing (face-centred cubic gridj^. In Figs. |4j&l5] we plot the re- 
sulting contours of ar/c^k~ over wave number \k\ and smoothing 
scale h (both normalised by the nearest-neighbour distance dnn) or 
Nh on the right axes (except for the Gaussian kernel when Nh is ill- 
defined and we give N, instead) for two wave directions, one being 
a nearest-neighbour direction. 



3.2. 1 Linear stability against pairing 

The top sub-panels of Figs.|4l&[5]refer to the longitudinal eigenvalue 
co^, when green and red contours are for, respectively, oj^^ x c^k^ and 
ai^ < 0, the latter indicative of the pairing instability. For the Gaus- 
sian kernel (truncated at I6it; Fig.|4l( a>f^ >0 eveiywhere, proving 
its stabilit)0, similar to the H0CT4 and, in particular the higher- 
degree, Wendland kernels. In contrast, all the B-spline kernels ob- 
tain < at sufficiently large A'^. 

The quintic spline, Wendland C^, and HOCT4 kernel each 
have a region of (^i^ < for \k\ close to the Nyquist frequency and 
Nh ~ 100, Nh a 40, and A'^ ~ 150, respectively. In numerical exper- 
iments similar to those described in 34.11 the corresponding insta- 
bility for the quintic spline and Wendland C" kernels can be trig- 
gered by very small random perturbations to the grid equilibrium. 



' Avoiding an obviously unstable configuration, such as a cubic lattice, 
which may result in (t)^ < simply because the configuration itself was un- 
stable, not the numerical scheme. 

There is in fact o)^ < at values for /; larger than plotted. In agreement 
with our analysis in j|3.1.2l this is caused by truncating the Gaussian, which 
(like any other modification to avoid infinite neighbour numbers) invali- 
dates the non-negativity of its Fourier transform. These theoretical results 
are confirmed by numerical findings of D. Price (referee report), who re- 
ports pairing at lai'ge /i/rf„„ for the truncated Gaussian. 
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Figure 6. Ratio of SPH sound speed csph = for the equation of state 

P ocp^/^ to the correct c for the kemel-Nn combinations of Table|2]and three 
difterent wave directions (as indicated) for paiticles in densest-sphere pack- 
ing with nearest-neighbour separation These curves are horizontal cuts 
through the top panels of Figs.|4]&|5] The thin vertical lines indicate sound 
with wavelength A = 8/i, corresponding to the dashed lines in Figs.|4]s|5] 



However, such modes are absent in glass-like configurations, which 
naturally emerge by 'cooling' initially random distributions. This 
strongly suggests, that these kemel-Nn combinations can be safely 
used in practice. Whether this also applies to the H0CT4 kernel at 
Nh ~ 150 we cannot say, as we have not run test simulations for this 
kernel. Note, that these islands of linear instability at small Nh are 
not in contradiction to the relation between kernel Fourier trans- 
form and stability and are quite different from the situation for the 
B-splines, which are only stable for sufficiently small Nh. 



3.2.2 Linear transverse instability? 

The bottom sub-panels of Figs. |4]S[5] show cj^^lc^k^, when both 
families of kernels have ~ Iw^,,! ^ c^k' with either sign oc- 
curring. < implies growing transverse modej^l which we 
indeed found in simulations starting from a slightly perturbed 
densest-sphere packing. However, such modes are not present in 
glass-like configurations, which strongly suggests, that transverse 
modes are not a problem in practice. 



" iRead et al] )2O10l) associate oj^j 2 ^ " v^'Ah a 'banding instability' which 
appeared near a contact discontinuity in some of their simulations. How- 
ever, they fail to provide convincing arguments for this connection, as their 
stability analysis is compromised by the use of the unstable cubic lattice. 
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kernel 


Nh 


Ni, 


' '/ "-'nil 


hibltn) 


cubic spline 


42 


6.90 


1.052 


1.181 


cubic spline 


55 


9.04 


1.151 


1.292 


quartic spline 


60 


7.29 


1.072 


1.203 


quintic spline 


180 


17.00 


1.421 


1.595 


Wendland C' 


100 


13.77 


1.325 


1.487 


Wendland C* 


200 


18.58 


1.464 


1.643 


Wendland 


400 


27.22 


1.662 


1.866 


HOCT4 


442 


42.10 


1.923 


2.158 


Gaussian 


oo 


10.00 


1.191 


1.337 


Gaussian 


DO 


20.00 


1.500 


1.684 



quartic spline, W^ = 200 



Wendland C^, N =200 



Table 2. Some quantities (defined in j)2.U for kernel-A'/f combinations used 
in Fig.|6]and the test simulations of ^ t/„n is the nearest-neighbour distance 
for densest-sphere packing, which has number density n = 'J2d^^ . The cu- 
bic spline with Nn x 42 is the most common choice in astrophysical sim- 
ulations, the other Nh values for the B-spline are near the pairing-stability 
boundary, hence obtaining close to the greatest possible reduction of the 
'Eo errors'. For Nh = 100, 200, 400, we picked the Wendland kernel which 
gave best results for the vortex test of j|4.2l 



3.3 Numerical resolution of sound waves 

The dashed lines in Figs.|4]&[5]indicate sound with wavelength A = 
8/7. For h > d„„, such sound waves are well resolved in the sense that 
the sound speed is accurate to < 1 %. This is similar to grid methods, 
which typically require about eight cells to resolve a wavelength. 

The effective SPH sound speed can be defined as Csph = t^w/lk]. 
In Fig.|6]we plot the ratio between csph and the coiTect sound speed 
as function of wave number for three different wave directions and 
the ten kernel- A'h combinations of Table [2] (which also gives their 
formal resolutions). The transition from csph ~ c for long waves to 
Csph ^ c for short waves occurs at \k\d,,„ > 1, but towards longer 
waves for larger h/dj^„, as expected. 

For resolved waves (A > Sh: left of the thin vertical lines in 
Fig-Ell, t'spH obtains a value close to c, but with clear differences be- 
tween the various kemel-Nn combinations. Surprisingly, the stan- 
dard cubic spline kernel, which is used almost exclusively in astro- 
physics, performs very poorly with errors of few percent, for both 
N/j = 42 and 55. This is in stark contrast to the quartic spline with 
similar Nu = 60 but Csph accurate to < 1%. Moreover, the quartic 
spline with Nu = 60 resolves shorter waves better than the cubic 
spline with a smaller = 55, in agreement with Table|2] 

We should note that these results for the numerical sound 
speed assume a perfectly smooth simulated flow. In practice, parti- 
cle disorder degrades the performance, in particular for smaller N^, 
and the resolution of SPH is limited by the need to suppress this 
degradation via increasing h (and A'^). 



4 TEST SIMULATIONS 

In order to assess the Wendland kernels and compare them to the 
standard B-spline kernels in practice, we present some test sim- 
ulations which emphasise the pairing, strong shear, and shocks. 
All these simulations are done in 3D using periodic boundary 
conditions, P = Kp^l^, conservative SPH (equation [Sj, and the 
ICullen & DehnerJ hoid) artificial viscosity treatment, which in- 
vokes dissipation only for compressive flow s , and an artificial con- 
ductivity similar to that o f iRead & Havfieldl ( 1201 2h . For some tests 
we used various values of Nh per kernel, but mostly those listed in 
Tabled 
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Figure 7. Final .t and y positions for particles at |z| < c/nn /2 for the tests of 
34.11 Symbol size is linear in z'- two overlapping symbols of the same size 
indicate particle pairing. 



4.1 Pairing in practice 

In order to test our theoretical predictions regarding the pairing in- 
stability, we evolve noisy initial conditions with 32000 particles 
until equilibrium is reached. Initially, x, = 0, while the initial jc, 
are generated from densest-sphere packing by adding normally dis- 
tributed offsets with (ID) standard deviation of one unperturbed 
nearest-neighbour distance d„„. To enable a uniform-density equi- 
librium (a glass), we suppress viscous heating. 

The typical outcome of these simulations is either a glass-like 
configuration (right panel of Fig. [7) or a distribution with particle 
pairs (left panel of Fig. [7}. In order to quantify these outcomes, we 
compute for each particle the ratio 



--mm[\Xi-Xj\]IHi 

j*' 



(30) 



between its actual nearest-neighbour distance and kernel- support 
radius. The maximum possible value for r^jn occurs for densest- 
sphere packing, when \x; -Xj\^ dnn = («/V2)''^ with n the number 
density. Replacing p, in equation ^ with m,7i, we obtain 



,<{3NH/2'"ny'\ 



Thus, the ratio 



_ ^ mill,/ . J -^/l 

- (3yV„/25/2;r)>/3 " "/iP 1 -d~; 



(31) 



(32) 



is an indicator for the regularity of the particle distribution around 
particle i. It obtains a value very close to one for perfect densest- 
sphere packing and near zero for pairing, while a glass typically 



gives q„ 



,~0.7. 



Fig. [8] plots the final value for the overall minimum of gniin.i 
for each of a set of simulations. For all values tested for Nh (up 
to 700), the Wendland kernels show no indication of a single par- 
ticle pair. This is in stark contrast to the B-spline kernels, all of 
which suffer from particle pairing. The pairing occurs at Nh > 67 
and 190 for the quartic, and quintic spline, respectively, whereas for 
the cubic spline mmj{q^i„ j] approaches zero more gradually, with 
min,{i5',^i„ , j < 0.16 at A';; > 55. These thresholds match quite well 
the suggestions of the linear stability analysis in Figs. |4lS|5] (except 
that the indications of instability of the quintic spline at Nh ~ 100 
and the Wendland kernel at Nh ~ 40 are not reflected in our tests 
here). The quintic (and higher-order) splines are the only option 
amongst the B-spline kernels for A';; appreciably larger than ~50. 

We also note that min/jgnjin,) grows substantially faster, in par- 
ticularly early on, for the Wendland kernels than for the B-splines, 
especially when operating close to the stability boundary. 



© 2012 RAS, MNRAS OQO.fnfTS] 



10 Walter Dehnen & Hossam Aly 




0.1 0.2 „ 0.3 0,4 0.1 0.2 „ 0.3 0,4 0.1 0.2 „ 0.3 0.4 0.1 0.2 „ 0.3 0.4 

ri r\ ri ri 

Figure 9. Velocity profiles the Gresho-Chan vortex test at time t=\ with the lowest resolution of Wid = 50.8 for three different kemel-A'H pairs (as indicated) 

and one simulation {right panel) with Pq = in equation (33a), when the 'Eq en'or' is naturally much weaker. Only particles at \z\ < 0.05 are plotted. 
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Figure 8. Final value of the overall minimum of q^i^j (equation 1321 for 
simulations starting from noisy initial conditions for all kernels of Table [T] 
(same colour coding as in Fig. [3] using the density coiTection of j|2.5l for the 
Wendland kernels) as function of Nu- For densest-sphere packing 5„,i„ , » 1, 
for a glass <s: ijminj < 1, while ^minj ~ indicates particle pairing. 



4.2 The Gresho-Chan vortex test 

As discussed in the introduction, particle disorder is unavoidably 
generated in shearing flows, inducing 'Eq errors' in the forces and 
causing modelling errors. A critical test of this situation consists 
of a difl' erentially rotating fluid of unifo rm density in centrifugal 
balance jCresho & Chanl ll990l. see also iLiska & Wendroflj |2003| . 



ISpringell^ld and lRead & Havfieldl2012ir The pressure and az- 
imuthal velocity are 



P = 



Po + l2.5R^+4-20R + 4\n(5R) 
[/'o + 2(21n2-l) 

(5R forO <R<0.2, 
^ 2-5R for 0.2 < < 0.4, 
for 0.4 < R 



for < < 0.2, 

for 0.2^ R< 0.4, (33a) 

for OA^R; 

(33b) 



with Po = 5 and R the cylindrical radius. We start our simulations 
from densest-sphere packing with effective one-dimensional parti- 
cle numbers A'id = 51, 102, 203, or 406. The initial velocities and 
pressure are set as in equations (133) . 

There are three different causes for errors in this test. First, an 
ov erly viscous me thod reduces the differential rotation, as shown 
bv lSpringellj2010l) : this effect is absent from our simulations ow- 
ing to the usage of the lCullen & DehnenI ( |2010|) dissipation switch. 
Second, the 'Eq eiTor' generates noise in the velocities which in 
turn triggers some viscosity. Finally, finite resolution implies that 
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Figure 10. Convergence of the L[ velocity eiTor for the Gresho-Chan test 
with increasing number of particles (for a cubic lattice Nid equals the num- 
ber of cells along one side of the computational domain) for various kemel- 
N[j combinations (those with filled squares are shown in Fig.|9]- A compari- 
son with Fig. 6 of Read & Hayfield (20 1 2) shows that the Wendland ker- 
nel with Nh = 400 performs better than the H0CT4 kernel with Nh = 442. 



the sharp velocity kinks at ^ = 0.2 and 0.4 cannot be fully resolved 
(in fact, the initial conditions are not in SPH equilibrium because 
the pressure gradient at these points is smoothed such that the SPH 
acceleration is not exactly balanced with the centrifugal force). 

In Fig.|9]we plot the azimuthal velocity at time f = 1 for a sub- 
set of all particles at our lowest resolution of A'lo = 51 for four dif- 
ferent kemel-A'// combinations. The leftmost is the standard cubic 
spline with = 55, which considerably suffers from particle dis- 
order and hence 'Eo errors (but also obtains too low atR< 0.2). 

The second is the Wendland C* kernel with Nh = 100, which 
still suffers from the 'Eo en'or'. The last two are for the Wendland 
C'' kernel with Nh = 400 and the Wendland kernel with Nh = WO 
but with /'o = in equation l |33at . In both cases, the 'Eq error' is 
much reduced (and the accuracy limited by resolution) either be- 
cause of large neighbour number or because of a reduced pressure. 

In Fig. [Tol we plot the convergence of the Li velocity error 
with increasing numerical resolution for all the kernels of Table [T] 
but with another Nh for each, see also Table|2] For the B-splines, we 
pick a large Nh which still gives sufficient stability against pairing, 
while for Nh = 100, 200, and 400 we show the Wendland kernel 
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cubic spline N„ = 55 quartic spline JV„ = 60 Wendland JV„= 100 quintic spline 180 Wendland C JV„=200 Wendland C JV„=400 
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Figure 11. Velocity, density, and thermal energy profiles for the shock tube test at ^ = 0.2 for the same kemel-Nu combinations as in Fig. llOl in order of 
increasing Nu (points: particles at lyl, |z| < 0.005) and the exact solutions (solid). The Li velocity en'ors reported are computed for the range -0.4 <x< 0.5. 



that gave best results. For the cu bic spUne, the re sults agree with 
the no-viscosity case in Fig. 6 of ISpringel ( l2010l) . demonstrating 
that our dissipation switch effectively yields inviscid SPH. We also 
see that the rate of convergence (the slope of the various curves) is 
lower for the cubic spline than any other kernel. This is caused by 
systematically too low in the rigidly rotating part at R< 0.2 (see 
leftmost panel if Fig.[9]( at all resolutions. The good performance of 
the quartic spline is quite sutprising, in particular given the rather 
low Nh- The quintic spline at A'h = 180 and the Wendland kernel 
at Nh = 200 obtain very similar convergence, but are clearly topped 
by the Wendland kernel at Njj = 400, demonstrating that high 
neighbour number is really helpful in strong shear flows. 



4.3 Shocks 



Our final test is the classical lsodi ( Il978l) shock tube, a ID Riemann 
problem, corresponding to an initial discontinuity in density and 
pressure. Unlike most published applications of this test, we per- 
form 3D simulations with glass-like initial conditions. Our objec- 
tive here is (1) to verify the 'Eo-error' reductions at larger Nu and 
(2) the resulting trade-off with the resolution across the shock and 
contact discontinuities. Other than for the vortex tests of 34.21 we 
only consider one value for the number A' of particles but the same 
six kemel-A'/y combinations as in Fig. (TO) The resulting profiles of 
velocity, density, and thermal energy are plotted in Fig. Ill [together 
with the exact solutions. 

Note that the usual over-shooting of the thermal energy near 
the contact discontinuity (at x = 0.17) is prevented by our artifi- 
cial conductivity treatment. This is not optimised and likely over- 
smoothes the thermal energy (and with it the density). However, 
here we concentrate on the velocity. 

For the cubic spline with Nu = 55, there is significant ve- 
locity noise in the post-shock region. This is caused by the re- 
ordering of the particle positions after the particle distribution be- 
comes anisotropically compressed in the shock. This type of noise 
is a well-kn own issue with multi-dimen sional SPH simulations of 
shocks (e.g. ISpringelllioTol : |Pricell2012l) . With increasing Nh the 
velocity noise is reduced, but because of the smoothing of the ve- 
locity jump at the shock (at x = 0.378) the Li velocity error does 
not approach zero for large Nh- 



Instead, for sufficiently large Nh (> 200 in this test), the Li 
velocity error saturates: any 'Eo-error' reduction for larger Nh is 
balanced by a loss of resolution. The only disadvantage of larger Nh 
is an increased computational cost (by a factor ~ 1.5 when moving 
from the quintic spline with Nh = 180 to the Wendland kernel 
with A'^ = 400, see Fig.fTIIl. 



5 DISCUSSION 

5.1 What causes the pairing instabiUty? 

The Wendland kernels have an inflection point and yet show no 
signs of the pairing instability. This clearly demonst rates that the 
tradit ional ideas for the origin of this instability (a la lSwegle et al.l 
1 19951 see the introduction) were incorrect. Instead, our linear sta- 
bility analysis shows that in the limit of large Nh pairing is caused 
by a negative kernel Fourier transform W, whereas the related ten- 
sile instability with the same symptoms is caused by an (effective) 
negative pressure. While it is intuitively clear that negative pres- 
sure causes pairing, the effect of < is less obvious. Therefore, 
we now provide another explanation, not restricted to large Nh. 



5.1.1 The pairing instability as artifact of the density estimator 

By their derivation from the Lagrangian (|2), the SPH forces ;7!,x, = 
dJl/dXj = -dU/dXj tend to reduce the estimated total thermal en- 
ergy U = Yji'^if(Pi,-'^i) at fixed entrop\F^ Thus, hydrostatic equi- 
librium corresponds to an extremum of U, and stable equilibrium 
to a minimum when small positional changes meet opposing for- 
ces. Minimal IJ is obtained for uniform p,, since a re-distribution 
of the particles in the same volume but with a spread of p, gives 
larger U (assuming uniform 5,). An equilibrium is meta-stable, if 
U is only a local (but not the global) minimum. Several extrema 
can occur if different particle distributions, each obtaining (near-) 
uniform p,, have different average p. Consider, for example, parti- 
cles in densest-sphere packing, replace each by a pair and increase 



'~ This holds, of course, also for an isothermal gas, when u is constant, but 
not the entropy .v, so that (du/dp)s 0. 
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the spacing by 2''^, so that the average density p (but not p) remains 
unchanged. This fully paired distribution is in equilibrium with uni- 
form Pi, but the effective neighbour number is reduced by a factor 2 
(for the same smoothing scale). Now if p(Nh/2) < p(Nh), the paired 
distribution has lower Lf than the original and is favoured. 

In practice (and in our simulations in 34. U . the pairing insta- 
bility appears gradually: for Nu just beyond the stability boundary, 
only few particle pairs form and the effective reduction of Nh is by 
a factor / < 2. We conclude, therefore, that 

pairing occurs if p{Nul f) < piNu) for some 1 < / ^ 2. 

From Fig.[3]we see that for the B-spline kernels piNn) always has a 
minimum and hence satisfies our condition, while this never occurs 
for the Wendland or H0CT4 kemelsB The stability boundary (be- 
tween squares and crosses in Fig. |3) is towards slightly larger Nh 
than the minimum of p(Nu), indicating / > 1 (but also note that the 
curves are based on a regular grid instead of a glass as the squares). 

5.1.2 The relation to 'Eg errors' and particle re-ordering 

A disordered particle distribution is typically not in equilibrium, but 
has non-uniform p; and hence non-minimal U. The SPH forces, in 
particular their 'Eo errors' (which occur even for constant pressure), 
then drive the evolution towards smaller U and h ence equilib rium 
with either a glass-like order or pairing (see also |Pricell2012l . §5). 
Thus, the minimisation of U is the underlying driver for both the 
particle re-ordering capability of SPH and the pairing instability. 
This also means that when operating near the stability boundary, 
for example using A'^ = 55 for the cubic spline, this re-ordering is 
much reduced. This is why in Fig. [8] the transition between glass 
and pairing is not abrupt: for Nh just below the stability boundary 
the glass-formation, which relies on the re-ordering mechanism, is 
very slow and not finished by the end of our test simulations. 

An immediate coroUaiy of these considerations is that any 
SPH-like method without 'Eq errors' does not have an automatic 
re-ordering mechanism. This applies to modifications of the force 
equation that a v oid th e 'Eo error', but also to the method of 
iHeB & Springeil tOldt . which employs a Voronoi tessellation to 
obtain the density estimates p, used in the particle Lagrangian (|2)- 
The tessellation constructs a partition of unity, such that different 
particle distributions with uniform p, have exactly the same average 
p, i.e. the global minimum of U is highly degenerate. This method 
has neither a pairing instability, nor 'Eo errors', nor the re-ordering 
capacity of SPH, but requires additional terms for that latter pur- 
pose. 

5.2 Are there more useful kernel functions? 

Neither the B-splines nor the Wendland functions have been de- 
signed with SPH or the task of density estimation in mind, but de- 
rive from interpolation of function values y, for given points x, . 

The B-splines were constructed to exactly interpolate polyno- 
mials on a regular ID grid. However, this for itself is not a desirable 
property in the context of SPH, in particular for 2D and 3D. 

The Wendland functions were designed for interpolation of 
scattered multi-dimensional data, viz 

six) = 'Zja;jw(\x-Xj\). 

The density correction of 32.5l does not affect these arguments, because 
during a simulation e in equation is fixed and in terms of our consider- 
ations here the solid curves in Fig.|3]are simply lowered by a constant. 




Figure 12. Wall-clock timings for an (average) single SPH time step (using 
four processors) for N = particles as function of neighbour number Nh 
for the kernels of Table [T] 

The coefficients Oj are determined by matching the interpolant s{x) 
to the function values, resulting in the linear equations 

yi = lLjajV;<\Xi-x-^, i= l...n. 

If the matrix W,^ = w(\Xi - Xj\) is positive definite for any choice 
of n points Xj, then this equation can always be solved. Moreover, 
if the function w{r) has compact support, then W is sparse, which 
greatly reduces the complexity of the problem. The Wendland func- 
tions were designed to fit this bill. As a side effect they have non- 
negative Fourier transform (according to lBochnetlll933l) , which to- 
gether with their compact support, smoothness, and computational 
simplicity makes them ideal for SPH with large Nh. 

So far, the Wendland functions are the only kernels which are 
stable against pairing for all Nh and satisfy all other desirable prop- 
erties from the list on page|3] 

5.3 What is the SPH resolution scale? 

In smooth flows, i.e. in the absence of particle disorder, the only 
error of the SPH estimates is the bias induced by the smoothing 
operation. For example, assuming a smooth density field 

p, »p(x,) + \cP- V'-pixd+OQi") (34) 

(e.g. lMonagha3ll985l :[ Silverma3l986l) with cr defined in equation 
([8}. Since cr also sets the resolution of sound waves ( 33.1.3K our 
definition dlOb . h = 2a, of the SPH resolution scale is appropriate 
for smooth flows. The result I l34t is the basis for the traditional 
claim of (9(/i^) convergence for smooth flows. True flow disconti- 
nuities are smeared out over a length scale comparable to h (though 
we have not made a detailed investigation of this). 

In practice, however, particle disorder affects the performance 
and, as our test simulations demonstrated, the actual resolution of 
SPH can be much worse than the smooth-flow limit suggests. 

5.4 Are large neighbour numbers sensible? 

There is no consensus about the best neighbour number in SP H: tra- 
dition ally the cubic spline kernel is used with Nh ~ 42, while |Pric3 
( I2OI2I) fav ours Nh = 57 (at o r even beyond the pairing-instability 
limit) and iRead et al. I ( I2OIOI) use their HOCT4 kenrel with even 
Nh = 442 (corresponding to a 1.7 times larger h). From a pragmatic 
point of view, the number N of particles, the neighbour number Nh, 
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Figure 13. As Fig. 1101 except that the ,v-axis shows the computational costs. 

and the smoothing kernel (and between them the numerical resolu- 
tion) are numerical parameters which can be chosen to optimise 
the efficiency of the simulation. The critical question therefore is: 

Which combination of A' and Nh (and kernel) most ef- 
ficiently models a given problem at a desired fidelity? 

Clearly, this will depend on the problem at hand as well as the de- 
sired fidelity. However, if the problem contains any chaotic or tur- 
bulent flows, as is common in star- and galaxy formation, then the 
situation exemplified in the Gresho-Chan vortex test of 94.2l is not 
atypical and large may be required for sufficient accuracy. 

But are high neighbour numbers affordable? In Fig. [T2] we 
plot the computational cost versus Nh for different kernels. At 
Nh S 400 the costs rise sub-linearly with Nh (because at low Nh 
SPH is data- rather than computation-dominated) and high Nh are 
well affordable. In the case of the vortex test, they are absolutely 
necessary as Fig.[T3]demonstrates: for a given numerical accuracy, 
our highest A';^ makes optimal use of the computational resources 
(in our code memoiy usage does not significantly depend on Nh, so 
CPU time is the only relevant resource). 



6 SUMMARY 

Particle disorder is unavoidable in strong shear (ubiquitous in astro- 
physical flows) and causes random errors of the SPH force estima- 
tor. The good news is that particle disorder is less severe than Pois- 
sonian shot noise and t he resulting force errors (which are domi- 
nated by the 'Eo' term of lRead et alj201(il) are not catastrophic. The 
bad news, however, is that these errors are still significant enough 
to spoil the convergence of SPH. 

In this study we investigated the option to reduce the 'Eo er- 
rors' by increasing the neighbour number in conjunction with a 
change of the smoothing kernel. Switching from the cubic to the 
quintic spline at fixed resolution h increases the neighbour number 
Nh only by a factoi0 1.74, hardly enough to combat 'Eo errors'. 



Using our definition l IlOl for the smoothing scale h. The conventional 
factor is 3.375, almost twice 1.74, but fomially effects to a loss of resolution, 
since the conventional value for h of the B-spline kernels is inappropriate. 



For a significant reduction of the these errors one has to trade reso- 
lution and significantly increase Nh beyond conventional values. 

The main obstacle with this approach is the pairing instability, 
which occurs for large Nh with the traditional SPH smoothing ker- 
nels. In ^and appendix [A] we have performed (it appears for the 
first time) a complete linear stability analysis for conservative SPH 
in any number of spatial dimensions. This analysis shows that SPH 
smoothing kernels whose Fourier transform is negative for some 
wave vector k will inevitably trigger the SPH pairing instability 
at sufficiently large neighbour number Nh- Such kernels therefore 
require Nh to not exceed a certain threshold in order to avoid the 
pairing instability (not to be confused with the tensile instability, 
which has the same symptoms but is caused by a negative effective 
pressure independent of the kernel properties). 

Intuitively, the pairing instability can be understood in terms 
of the SPH density estimator: if a paired particle distribution ob- 
tains a lower average estimated density, its estimated total thermal 
energy U is smaller and hence favourable. Otherwise, the smallest 
U occurs for a regular distribution, driving the automatic mainte- 
nance of particle order, a fundamental ingredient of SPH. 

The IWendlan^ ( Il995h functions, presented in 32.2.31 have 
been constructed, albeit for different reasons, to possess a non- 
negative Fourier transform, and be of compact support with simple 
functional form. The first property and the findings from our tests 
in 34. II demonstrate the remarkable fact that these kernels are sta- 
ble against pairing for all neighbour numbers (this disproves the 
long-cultivated myth that the pairing instability was caused by a 
maximum in the kernel gradient). Our 3D test simulations show 
that the cubic, quartic, and quintic spline kernels become unstable 
to pairing for Nh > 55, 67, and 190, respectively (see Fig. [8), but 
operating close to these thresholds cannot be recommended. 

A drawback of the Wendland kernels is a comparably large 
density en'or at low Nh- As we argue in 35.1.11 this error is directly 
related to the stability against pairing. However, in 32.5l we present 
a simple method to correct for this error without affecting the sta- 
bility properties and without any other adverse effects. 

We conclude, therefore, that the Wendland functions are ideal 
candidates for SPH smoothing kernels, in particular when large Nh 
are desired, since they are computationally superior to the high- 
order B-splines. All other alternative kernels proposed in the liter- 
ature are computationally more d emanding and are either centrally 
spiked, like the H0CT4 kernel o flRead et aljj201(]|). or susceptible 
to pairing like the B-splines (e.g. ICabezon et al.ll2008l) . 

Our tests of Section |4] show that simulations of both strong 
shear flows and shocks benefit from large Nh. These tests suggest 
that for Nh ~ 200 and 400, respectively, the Wendland and 
kernels are most suitable. Compared to Nh = 55 with the standard 
cubic spline kernel, these kemel-A'/y combinations have lower res- 
olution (h increased by factors 1.27 and 1.44, respectively), but ob- 
tain much better convergence in our tests. 

For small neighbour numbers, however, these tests and our 
linear stability analysis unexpectedly show that the quartic B- 
spline kernel with Nh = 60 is clearly superior to the traditional 
cubic spline and can compete with the Wendland kernel with 
Nh = 100. The reason for this astonishing performance of the quar- 
tic spline is unclear, perhaps the fact that near x = this spline is 
more than three times continuously differentiable plays a role. 

We note that, while the higher-degree Wendland functions 
are new to SPH, the Wendland kernel has already been used 
jMonaghai] 1201 ll . for example, employs it for 2D simulations). 
However, while its immunity to the pairing instability has been 
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noted (e.g. iRobinsonI |200^ we are not aware of any explana- 
tion (previous to ours) nor of any other systematic investigation of 
the suitabihty of the Wendland functions for SPH. 
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APPENDIX A: LINEAR STABILITY ANALYSIS 

We start from an equilibrium with particles of equal mass m on a 
regular grid and impose a plane-wave perturbation to the unper- 
turbed positions i, (a bar denotes a quantity obtained for the unper- 
turbed equilibrium): 

x,.=je,+f,., f,. = aO,(jt:,f), <!),(■«, f) = e''*'*'"'"'. (Al) 

as in equation l |20t . We derive the dispersion relation (jj(k,a) by 
equating the SPH force imposed by the perturbation (to linear or- 
der) to its acceleration 

l = -ao?'b,{x,t). (A2) 

To obtain the perturbed SPH forces to linear order, we develop the 
internal energy of the system, and hence the SPH density estimate, 
to second order in f,. If p, =p+p,\ +pa with p,i and pa the first and 
second-order density corrections, respectively, then 

(A3) 

AI Fixed h 

Let us first consider the simple case of constant h = h which remains 
unchanged during the perturbation. Then, 

p,„ =£.,„/«! with Q,„ = Y,jm(^,yV)"W(x,jJi). (A4) 
Inserting 

^,,=^,-^-a(l-e-'"'n(l), (A5) 
into JA4b gives 

fti =0, Y.iin{l-e'''"^'i)a-VW{Xij,K) = i^ia-t, (A6) 
where (assuming a symmetric particle distribution) 
t{k) = Y,jm?.mk-XjVW(Xj,h)- (A7) 
We can then derive 

'LtQkiden/d^i = T^teki I,jmW/(xi,j,h){Sik-6ij) 
= Zj'"{gn+gji)'^W{Xij,h) 
= i<S>ia t ( 1 + e"'*-*") ^W{Xij, h) 
= <t>iatt, (A8) 
i JlkdQn/d^, = j:,jm(^,yVWW(x,jM6it-6ij) 
= 2j:jm(^yV)VW{x,jJi) 
= 20, j:jm{\-e-"'-^")a-VVW{xtj,h) 
= 2<t>iaT (A9) 
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with 

T(A:) = 2; m ( 1 - cos k-Xj) V*-' W{Xj, h). 
Inserting these results into l |A3t , we get 



P P 



2P\ a tt 
— — 

P I P^ 



(AlO) 



(All) 



Inserting these into ( IA16t gives 



\2 + v I 



kk+- 



vk^a] 
y+2 



(A20) 



A2 Adaptive smoothing 

If the hj are adapted such that M, = h^pj remains a global constant 
M, the estimated density is simply /3, = Mhj". We start by expand- 
ing Mj to second order in both a and = ln(hj/h). Using a prime to 
denote differentiation w.r.t. Inh, we have 

M, = M,,+l7"gn +r]{h'p)' + fh^ga+vih^Qn)' + ^r^^h^p)" (A12a) 
= Mi,+h" [eii+TjvpCl+^ga + rie]^ +T]vgn + ^rfv^P^] (A12b) 

with n 1 as defined in equation l|4) and 
1 d~{h"p) 1 dih'pQ) 



n 



v^h'p d(\nh)^ vh"p dlnh 



1. 



(A 13) 



By demanding M, = M, we obtain for the first- and second-order 
contributions to 77. 



en 

-vprin = H = — H r h r 

man 2q 



en (en)' 



eh__^_en^ 



(A14a) 
(A14b) 
(A 14c) 



2Q. 2vpfi2 pQ2 2pn3 

From these expressions and p; = pe""'' we obtain the first and sec- 
ond order density corrections 

n ' 

ea 



Pil = -ypin 



Pa=-^P1a + \y^prin 



(elr ^em-ii) 



(A15a) 
(A15b) 



2n 2vpQ? 2pQ? 

Inserting these expressions into equation ( IA3t we find with rela 
tions ( |A8t and (|A5} 

1 d{a-tt) 



P al - a tt 

£: = 2^+H^^ 

P\ pQ. p2n2 yp2Ql dXnh 

I , 2P\ a tt 

\ pIp^o? 



where 

^ n 



^_ i ain(pQ) ^^ 



(A16) 



(A17) 



A3 The limit A: -» 

From equation ((S} /i'^-VWCjc,/!) = -a(//W')/ainA, such that 

- vp.fi, = i;^.m,x,yVH/(x,y,/i,), (A18a) 

v-p,n, = I],'«,(^,-V)^iy(x,y,/;,). (A18b) 

Assuming local spatial isotropy we then get in the limit — > 

t -Q.pk, (A19a) 

t' -v{n-Cl)pk, (A19b) 

2 1. (A19c) 



(2f2 + vn)p p) v(n-f2)p, 



2 + v 2(v-h2) 
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